Raa of J /if) near mid-rapidity in heavy ion collisions at ^/snn = 200 GeV 
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We build up a model to reproduce the experimentally measured Raa of J/ip near midrapidty in 
Au+Au collision at ^smn — 200 GeV. The model takes into account the J/ip suppression from the 
quark-gluon plasma and hadron gas as well as the nuclear absorption of primordial charmonia and the 
regeneration effects at the hadronization stage, and hence is a generalization of the two component 
model introduced by Grandchamp and Rapp. The improvements in this work are twofold; the 
addition of the initial local temperature profile and a consistent use of QCD NLO formula for 
both the dissociation cross section in the hadron gas and the thermal decay widths in the quark- 
gluon plasma for the charmonium states. The initial local temperature profile is determined from 
the assumption that the local entropy density is proportional to a formula involving the number 
densities of the number of participants and of the binary collisions that reproduces the multiplicities 
of charged particles at chemical freeze-out. The initial local temperature profile brings about a kink 
in the Raa curve due to the initial melting of J/ip. The initially formed fireball, composed of weakly 
interacting quarks and gluons with thermal masses that are extracted from lattice QCD, follows an 
isentropic expansion with cylindrical symmetry. The fit reproduces well the Au+Au as well as the 
Cu+Cu data. The same method is applied to predict the Raa expected from the Pb+Pb collision 
at LHC energy. 
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I. INTRODUCTION 

Ever since the seminal work by Matsui and Satz 
J/ijj has been investigated for a long time as a diagnos- 
tic tool to probe the properties of hot nuclear matter 
created in the early stages of a relativistic heavy ion col- 
lision. The original claim stated that J/ip cannot be 
formed in the quark-gluon plasma(QGP) due to color 
screening between the charm and anti-charm quark and 
that such mechanism will lead to the suppression of J / ip 
production in a heavy ion collision if the QGP is formed. 
However, the situation was found to be more involved as 
recent lattice QCD results showed that J/ip might sur- 
vive past the critical temperature (T c ) and dissolve only 
at a higher temperature (Td ~ 2T C )0, Q. If so, it is im- 
portant to know the detailed temperature dependencies 
of the properties of the J/ip above T c , as large thermal 
width for example might still lead to a very small sur- 
vival rate of the J/-0JJ-Q. Moreover, if the quark-gluon 
plasma is formed and the number of charm quarks are 
large, there is an additional production mechanism that 
has to be estimated and comes from the formation of 
J/ip through the regeneration of charm and anti-charm 
that are well described by statistical approaches A 
successful phcnomcnological model that is able to explain 
the J/ip signals coming from p+p to heavy ion collisions 
at high energy should consistently include all the above 
mentioned ingredients. 
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The commonly measured and calculated nuclear modi- 
fication factor Raa shows whether the J/i p is suppressed 
or enhanced in the heavy ion collision [111 ]. For instance, 
if Raa = 1- the number of produced J/ip in A+A col- 
lision is equal to the number of produced J/ip in p+p 
collision at the same energy, times the number of bi- 
nary collisions of nucleons from the two colliding nuclei. 
If Raa < 1(> 1), the number of J/ip dissolved in hot 
nuclear matter is larger (smaller) than that regenerated 
from the quark-gluon plasma at the hadronization point. 
Presently Raa of J/i p is below unity up to the maximum 
energy of RHIC. 

Phenomenological models have been developed to de- 
scribe Raa of J/ip from SPS to RHIC, and to predict 
the upcoming LHC heavy ion data. One of the success- 
ful models to explain the heavy ion data is the thermal 
model [HI]. Here the J/ip production follows the statis- 
tical description at the chemical freeze out point, and 
thus can be thought of as being regenerated from the 
quark-gluon plasma at the hadronization point (~ T c ). 
In the so called two-component model by Grandchamp 
and Rapp(GR) one additionally takes into account 
the fact that the J/ip survives past T c and dissociates 
only at a higher temperature Td so that some of J/tp 
produced at initial stage survives the high temperature 
phase and reaches the final stage; these together with 
those regenerated from the quark-gluon plasma at the 
hadronization point constitutes the two components of 
the observed J/ip [HI]. In a third model the continuous 
dissociation and regeneration of J/ip is calculated using 
the Boltzman transport equation from Td to T c evolved 
using the hydrodynamics equation [l6j . 

Crucial quantities in the above models are the T^'s for 
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the charmonia; these are the melting temperatures for 
the Deby screened Coulomb type of charmonium bound 
states above T c . While, initial calculations seemed to 
suggest that Td ~ 2T C for J/tp and lower for the excited 
states [H, H, H3], recent studies claim that these could be 
lower due to the large thermal width the states acquire 
13, EH ■ The current experimental error for Raa of J / ip 
as a function of the number of participants are still too 
large to draw any definite conclusions about these differ- 
ent temperatures. However, one finds interesting charac- 
teristics if one takes the central values of the experimental 
data for the curve of Au+Au collision at Jsnn = 200 
GeV; there seems to be two sudden drops [l2j. One oc- 
curs initially at a very small number of participants and 
the other somewhere between 170 and 210. Assuming 
that the initial temperatures of the fireball is correlated 
to the number of participants, the two drops might be re- 
lated to two distinct initial temperatures. Indeed, as we 
will show in our two component model, the first sudden 
drop at small number of participants can be related to 
the dissociation of excited charmonia such as \c and ip'; 
the initial dissociation of these states will suppress the 
expected 30 to 40% feedback production of J/ ip from its 
excited states. The second sudden drop is found to be 
related to the initial dissociation of J/tp. In geometrical 
terms, the first drop occurs when there appears regions 
where the initial temperature is larger than the Td of the 
excited charmonia, and the second drop when it is larger 
than Td of J/tp. Moreover, if the thermal decay width of 
J ftp is not so large, these effects become prominent. 

In this work, we generalize the two-component model 
to take into account the initial hot region mentioned 
above, and calculate the Raa of J/ip at midrapidity in 
Au+Au collision at ^/snn = 200 GeV. We estimate the 
initial local temperature with the assumption that local 
entropy density is simply the linear combination of num- 
ber density of the number of participants and of binary 
collision, a formula that well reproduce the multiplicities 
of charged particles at chemical freeze-out stage in the 
scenario of isentropically expanding fireball. The aver- 
age temperature of the fireball is obtained from equating 
the entropy density to that of the quark-gluon plasma 
calculated using a non interacting gas of quarks and glu- 
ons with effective thermal masses extracted from lattice 
QCD calculation for the energy density and pressure [lj|. 
The baryon chemical potential is obtained from the ra- 
tio of the proton and the antiproton, and other chemical 
potentials from the conservation of respective flavors. 

In the two-component model, the primordial char- 
monia, which have not evolved to the final charmonia 
yet, first undergo nuclear absorption. Then after they 
form into the initially produced charmonia, they undergo 
through thermal decay in the quark-gluon plasma, and 
then hadronic decay in the hadron gas until the chem- 
ical freeze-out stage. On the other hand, the regener- 
ated charmonia only go through the hadronic decay in 
the hadron gas. Consequently, another important quan- 
tity in the model is the thermal widths of charmonia in 



the quark-gluon plasma and in hadron gas. We improve 
previous calculations by using the results obtained us- 
ing perturbative QCD up to the next to leading order 
(NLO) in the coupling constant [30| . In our work, this 
coupling constant is not a free parameter, because it is 
related to the screening mass of quark-gluon plasma and 
sequential melting temperatures of charmonia. The de- 
tails is mentioned in chapter VII. The relaxation factor 
of charm quarks in quark-gluon plasma is calculated up 
to leading order in perturbative QCD; this factor deter- 
mines the fraction of thermalized charm quarks and used 
to estimate the regeneration of charm quarks. 

The paper is organized as follows. In section II, we 
discuss nuclear absorption and introduce the necessary 
concepts. In section III, we introduce the initial local 
temperature profile. In section IV, we show how the ther- 
modynamic parameters are determined. In section V and 
VI we respectively discuss the thermal decay and regen- 
eration of charmonia. In the last two sections, we give 
the results and conclusions. In the appendix, we summa- 
rize the perturbative NLO formula for the dissociation 
cross section of charmonium. 



II. NUCLEAR ABSORPTION 

Within the Glauber model, heavy ion collision can be 
described with collisions of the nucleons. The model has 
two important scales, the number of participants and the 
number of binary collisions. Participants mean the nu- 
cleons in both colliding nuclei that go through inelastic 
scattering at least once. The rest are called spectators. 
Usually the amount of bulk matter created from heavy 
ion collision is proportional to the number of participants. 
The number of binary collisions counts only primary col- 
lisions of two nucleons. Usually hard particle such as 
the charmonium is produced through primary collision, 
because it requires large energy transfer; hence, its pro- 
duction number is proportional to the number of binary 
collisions. The two numbers are calculated in Glauber 
model as follows 12011: 
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where b is the impact parameter of the two colliding 
heavy nuclei, and A, B are their mass numbers. oy„ 
is the inelastic cross section of two nucleons, which is 
about 42 mb at ^/sWn = 200 GeV [ll[. T A ( B ) is called 
the thickness function defined as 



T 



A(B 



dzp A (B)(s, z), 



(2) 



3 



and Tab is called the overlap function. Here, Pa(b)(?) is 
the distribution function of one nucleon in nucleus A(B), 
for which we use the following Woods-Saxon model: 
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For the Au nucleus, r a = 6.38 fin, C = 0.535 fm [2l|]. The 
normalization factor p a is determined by the requirement 
/ d 3 rp(r) = 1. 

The charm and anticharm quark pair produced 
through primary collision is called the primordial char- 
monia before they form into on shell states. Such primor- 
dial charmonia can still be absorbed by nucleons passing 
through it. The survival rate at transverse position s 
from the nuclear absorption is 



dzdz' pa{s, z)p B (b — s, z') 



T AB (b) 

xexpl -(A-l) I dz A pA{s, z A )<J nuc 



xexp< - (B - 1) / dz B pB(b - s, z B )a nuc \. (4) 



Here, a nuc is the absorption cross section of primor- 
dial charmonia by a nucleon, and is about 1.5 mb at 
y/sNN = 200 GeV [HI . z and z' are the positions where 
the primordial charmonium is created in beam direction 
from the center of nucleus A and from the center of nu- 
cleus B respectively. The second line in Eq. ((U) is the 
absorption factor by nucleus A, and the third by nucleus 
B. Mass number A and B in the exponents of the sec- 
ond and third line are both subtracted by 1, because one 
nucleon of nucleus A and one nucleon of nucleus B take 
part in producing a charm pair, and cannot take part in 
absorbing the pair. 



III. INITIAL LOCAL TEMPERATURE 

To build up a model for the suppression, we will closely 
follow the space time picture of nucleus nucleus collision 
generally accepted from a hydrodynamic simulations (22j. 
After the heavy ions pass through each other, hot nu- 
clear matter is created with small net baryon density 
between the two receding nuclei. The hot matter is ex- 
pected to thcrmalize early; while the exact thermaliza- 
tion time could be controversial, we will just take the 
typical time accepted in hydrodynamic simulation given 
as r = 0.6fm/c[2lll|- 

The total entropy produced in heavy ion collision is es- 
timated from the multiplicity of produced particles. The 
multiplicity of charged particles near midrapidity was 
found to be well described under the assumption that 
it is proportional to a linear combination of the number 
of participants and the number of binary collisions. Be- 
cause both the multiplicity of particles and the entropy 
are extensive thermal quantities, it is assumed that the 



entropy per pseudorapidity is also proportional to the 
same linear combination of number of participants and 
binary collisions as the particle multiplicity. 
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In the linear combination for multiplicities of charged 
particles, x is 0.09 at ^/s^n = 130 GeV and 0.11 at 
\/s nn = 200 GeV; the same values are used for entropy 
[111 [26jj. The overall factor of 30.3 is determined so that 
multiplicities of charged particles are well reproduced un- 
der the condition that the ratio between entropy and mul- 
tiplicity is determined from the statistical model at the 
chemical freeze out point; c = S/M = s/iim, where S, M 
(s,nj\/) arc the total entropy and multiplicity (density) 
respectively. The upper figure of Fig. 2] shows that Eq. 
(|5"|). which is the solid line, reproduces the data well, when 
iscntropic expansion of hot nuclear matter is assumed. 

Dividing Eq. (JS|) by the volume, entropy density per 
pseudorapidity has the form 
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where n part , n co U ar e volume densities of the number 
of participants and binary collisions respectively. These 
densities at thermalization time r G , which is taken to be 
0.6 fm/c in this work, are given respectively as, 
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where s is the position vector on the transverse plane, and 
homogeneous densities in beam direction is assumed. 

Near midrapidity, the rapidity is approximated to the 
longitudinal velocity, 
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and is further assumed to be equal to the pseudorapidity. 
Supposing that the rapidity interval of interest is Ay, the 
total entropy within the interval is (dS/dr))Ay. Ignor- 
ing longitudinal acceleration, longitudinal size of nuclear 
matter within Ay is t (3 z w r Ay at thermalization time 
t q . Therefore, entropy density at thermalization time is 
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where A is the transverse area. 

The equation of state is required to extract the tem- 
perature from the entropy density. In the quasi-particle 
picture, the strongly interacting massless particles are 
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replaced by the noninteracting massive particles, whose 
properties are determined from the condition that they 
reproduce the thermal quantities such as energy density 
and pressure of the strongly interacting massless par- 
ticles. The thermal quantities of strongly interacting 
quarks and gluons are simulated by lattice QCD. In this 
work, we will use the parametrization for the effective 
thermal masses of quark and gluon extracted by Levai 
and Heinz from lattice QCD data for pure gauge, for 
Nt = 2 and for Nf = 4 [19| ■ They fix the degrees of free- 
dom of quarks and gluon to those in the vacuum, and use 
the forms of LO perturbative QCD for their masses, pa- 
rameterizing the the nonperturbative effects into a tem- 
perature dependent coupling g(T): 
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FIG. 1: Entropy density vs. temperature in the hadron gas 
(solid line) and in the quark-gluon plasma with parameters 
extracted from lattice QCD with two flavors (dotted line) and 
with four flavors (dashed line) 
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where 
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T c = 260, 140, 170 MeV and T c /A = 1.03, 1.03, 1.05 for 
pure gauge, Nf = 2 and Nf = 4 cases respectively. Here 
we choose the case of Nf = 4 for thermal masses, because 
its critical temperature is more reasonable than that of 
Nj = 2. 

With these effective thermal masses and the equation 
of state obtained by assuming a gas of weakly interact- 
ing quarks and gluons, the entropy density is related to 
temperature as given below, 
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The entropy density of quark-gluon plasma has a term 
proportional to the derivative of the squared thermal 
mass with respect to temperature, because thermal 
masses depend on temperature. But the term is to be 
canceled by the term induced from bag pressure to main- 
tain thermodynamic consistency (27j . The summation i 
includes gluon and quark of three flavors, although the 
lattice data is extracted for Nf — 4, we use the realistic 
number of flavors to get a more realistic magnitude for 
the entropy. The signs in the denominators are minus 
for gluon and plus for quarks. Dashed line and dotted 
line in Fig. [T] show the entropy density as a function of 
the temperature in quark-gluon plasma respectively for 
parameters extracted from lattice QCD for Nf = 4 and 
N f =2. 

From the entropy density of Eq. ((BJ) and the relation 
between entropy density and temperature given in Eq. 
(|13p. the initial local temperature distribution at r = t„ 




FIG. 2: Profiles of temperature along y-axis at various impact 
parameters 



is obtained as given in Fig. [2] and [3l y-axis is the direc- 
tion of impact parameter and x-axis is perpendicular to 
beam axis and y-axis. Fig. [5] shows temperature profiles 
along the y-axis at various impact parameters. As we can 
see, the maximum temperature becomes lower as impact 
parameter increases. Fig. [3] shows isothermal lines on 
transverse plane. In head-on collision case, the trans- 
verse radius of quark-gluon plasma is about 7.0 fm at 
thermalization time. This size decreases as the collision 
becomes peripheral. More importantly, the equithermal 
line of T = 380 MeV exists at b = 0, 4 fm, but does not 
at b = 8 fm. If the melting temperature of J/ip is 380 
MeV, the region where J/ip cannot be formed exists at 
b = and at b = 4 fm, but does not at b = 8 fm, which, 
as will be seen later, leads to a sudden drop to Raa curve 
of J ftp as & function of the number of participants. 
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FIG. 3: Isothermal lines on transverse plane at r = r with 
impact parameters, b = 0, 4, 8 fm 



IV. THERMAL QUANTITIES IN EXPANDING 
HOT NUCLEAR MATTER 



The hot nuclear matter created through a heavy ion 
collision expands in time; assuming iscntropic expansion, 
the entropy density and the temperature decreases along 
with it. Once the critical temperature is reached, the 
phase of the matter changes from quark-gluon plasma 
to the hadron gas. It is known that the phase transi- 
tion is a crossover at small baryon chemical potential 
and that there is a critical point where the phase transi- 
tion changes from a crossover to the first order as baryon 



chemical potential increases. While the exact location of 
the point is not certain yet, the phase transition is ex- 
pected to be still at the rapid crossover region at RHIC 
and at LHC energies. 

The statistical model is very successful in reproducing 
the observed particle ratios with two free parameters: the 
temperature and the baryon chemical potential at the 
chemical freeze-out stage. The other chemical potentials 
are obtained from the condition that each flavor must be 
conserved, which can be written as follows: 

^E n ^ = Z + N 
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Z -N 
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^E"A = o. 



(14) 



Here V is the volume of the hot nuclear matter of in- 
terest, Z and ./V the number of protons and neutrons in 
that volume respectively, and j the constituents of the 
matter: the constituents are the quarks and gluons in 
the quark-gluon plasma, while they are the mesons and 
baryons in the hadron gas. For hadron gas, all mesons 
and baryons whose masses are less than 1.5 GcV and 
2.0 GeV respectively are included in the sum. These 
upper limits for the hadron masses are taken such that 
inclusion of hadrons with masses higher than the upper 
bounds will not change the result; for example, if the 
upper limit of meson mass is set at 2.0 GeV, much more 
mesons are to be included, but their contribution is small. 
rij is the number density of particle j in the grandcanon- 
ical ensemble. Bj, I 3 j and Sj are baryon number, the 
third component of isospin and strangeness of particle j 
respectively. 

In our work, V will represent the volume in the midra- 
pidity region, which is defined by \y\ < 0.35. To deter- 
mine all the chemical potentials from Eq. (|14p . we have 
to know the two numbers Z and N. The two numbers 
are determined from two constraints: the first is that the 
ratio Z/N should be the same as in the original combined 
nuclei, and the second is that the observed ratio of pro- 
ton to antiproton within the rapidity is well reproduced. 
The second constraint needs some more explanation. In- 
stead of solving Eq. (TT4"|) for each N part , we will try to 
find the best fit for the experimentally observed p/p ratio 
as a function of N part by assuming that Z + N appearing 
in the right hand side of the first equation in Eq. (TT4")) is 
equal to the total number of participants in the collision 
divided by a constant number; the constant number is 
found to be 20.5. 

The relation between Z + N and the ratio of proton 
and antiproton is as following. If Z + N is very large, 
the nuclear matter with volume V will have much more 
baryons than antibaryons and thus baryon chemical po- 
tential will be large and positive; in this case, the proton 
to antiproton ratio will be much larger than 1. On the 
other hand, if Z + N is very small, the ratio will be close 
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FIG. 4: multiplicities of charged particles per pseudorapidity 
divided by a half number of participants (upper) and ratios of 
proton and antiproton (lower) near midrapidity at ^snn — 
200 GeV 



to 1. As one can see in the lower figure of Fig. [4] the 
best fit shown as the solid curve well reproduces the ex- 
perimental data. 

For the volume V, the following cylindrical form is used 
for simplicity: 



r + \ai_T 2 



(15) 



The factor 2 is multiplied to take into account the for- 
ward and the backward expansion in the beam direction. 
ft is the longitudinal velocity of nuclear matter, which 
is approximately equal to the rapidity near midrapidity. 
t is thermalization time, and r a is the initial radius of 
quark-gluon plasma on the transverse plane. As can be 
seen in Fig. [3J the shape of quark-gluon plasma is not 
a circle but almond-like except for the head-on collision 
case. 

In this work, the almond shape is transformed to a 
circle with the same area. a± is the transverse acceler- 
ation of quark-gluon plasma, which is set to 0.1 c 2 / fm 
[ToT ]. The acceleration is turned off, when the transverse 
velocity reaches 0.6 c. The initial total entropy of the 
quark-gluon plasma is calculated from a modified ver- 
sion of Eq. ([5]), where the N part and N co u arc not the 



total number of participants and total number of binary 
collisions but are respectively the number of participants 
and the number of binary collisions within quark-gluon 
plasma. Then, the entropy density is obtained by divid- 
ing the total entropy of the quark-gluon plasma by the 
initial volume as given in Eq. (j!5|) at r = tq. Finally, 
the entropy density at a later invariant time is obtained 
by assuming the expansion to be isentropic and that the 
volume expands as given in Eq. (|15[) , 

With the chemical potentials obtained from Eq. (|14l) . 
the entropy density determines temperature. The evo- 
lution of temperature of the quark-gluon plasma is ob- 
tained by equating the expression for the entropy density 
given in Eq. (|13[) to the previously determined entropy 
density at a later time. In contrast to the case for the 
quark-gluon plasma, it is assumed that the masses of all 
hadrons do not change at finite temperature. As we dis- 
cussed before, the order of the phase transition that we 
probe at RHIC seems to be a strong cross over. How- 
ever, the equation of states for the quark-gluon plasma 
and the hadron gas that we use gives a 1st order tran- 
sition. Nevertheless, we will still use these equation of 
states because they are simple, physically intuitive and 
easy to manipulate. Moreover, although at the physical 
quark masses, the order of phase transition is a rapid 
cross over, the change occurs at a very small temper- 
ature range so that effectively we can approximate the 
transition with a simple 1st order transition. 

The solid line in Fig. [TJ shows the correspondence 
between the entropy density to the temperature in the 
hadron gas. As can be seen in Fig. [TJ the entropy 
density of the hadron gas is higher than that of quark- 
gluon plasma at high temperature. These curves cross at 
around T = 210 MeV below which the entropy density 
of hadron gas is lower than that of quark-gluon plasma. 
But the curve of hadron gas and curves of quark-gluon 
plasma are very close to each other, which means that the 
phase transition from quark-gluon plasma to hadron gas 
is fast. As hot nuclear matter expands, its entropy den- 
sity and temperature decreases along the curve of quark- 
gluon plasma in Fig. [TJ At critical temperature, the 
nuclear matter transfers from the curve of quark-gluon 
plasma to that of hadron gas. If two curves meet at 
critical temperature and their derivatives with respect 
to temperature are the same, the phase transition is a 
crossover. If two curves meet but their derivatives are 
different, the phase transition is second order. But two 
curves in Fig. [TJ do not meet at critical temperature so 
the phase transition is first order. However, the phase 
transition takes short time, because two curves are close 
at critical temperature. This time interval is the period 
of mixed phase where two phases coexist. The upper 
figure of Fig. [S] shows time dependence of the average 
temperatures of the hot nuclear matter in head-on colli- 
sion. 

The lifetime of quark-gluon plasma phase in head-on 
Au+Au collision at ^/sjvw = 200 GeV is 5.63 fm/c. 
Mixed phase lasts for 1.44 fm/c. Temperature for chemi- 
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FIG. 5: The time evolution of average temperature(upper 
figure) and chemical potentials(lower figure) of hot nuclear 
matter at b = fm. The solid line, dotted line, and dashed 
line in the lower figure are baryon, strangeness and isospin 
chemical potentials respectively. 



cal freeze-out is set at 161 MeV, as motivated by thermal 
models. After chemical freeze-out, the number of parti- 
cles of each species is assumed to remain constant. Mul- 
tiplicities of charged particles and proton to antiproton 
ratio are all calculated at this temperature. Hadron gas 
phase lasts for 1.52 fm/c to chemical freeze-out. 

Charmonia produced outside hot region where ini- 
tial temperature is higher than dissociation tempera- 
ture of charmonia go through thermal decay caused by 
interactions with surrounding particles, which are the 
quarks and the gluons in the quark-gluon plasma, and the 
hadrons in the hadron gas or both in the mixed phase. 
Once thermal quantities of hot nuclear matter are given 
as functions of time, one needs to know the properties 
of charmonia in hot nuclear matter in order to calculate 
how many charmonia survive the thermal dissociation. 



r(T) 



d 3 k 



v rel (k) ni (k,T)af ss (k,T), (16) 



where i repents the quark or the gluon in the quark-gluon 
plasma, or the baryons and the meson in the hadron gas. 
di is the degeneracy factor of particle i, rii the number 
density of particle i in the grand canonical ensemble, v re i 
the relative velocity between the charmonium and the 
particle i, and af zss the dissociation cross section of char- 
monium by particle i. It is assumed that thermal width 
in the mixed phase is a linear combination of contribu- 
tions from the quark-gluon plasma and from the hadron 
gas as given in the following form: 



r = jt hg + (i - f)v^ GP , 



(17) 



where / is the fraction of the hadron gas in the mixed 
phase. 



A. dissociation cross section 

The crucial quantity in Eq. (|16p is the dissociation 
cross section. As emphasized before, one of the im- 
provements in our work over the previous two component 
model calculations is the consistent application of the 
NLO pcrturbative formula to calculate the dissociation 
of charmonia both in the quark-gluon plasma [2!| and in 
the hadron gas [3(j. In fact, it is difficult to describe the 
dissociation of charmonia both in the quark-gluon plasma 
and in the hadron gas with the same approach. As an ex- 
ample, GR in ref. 14| uses a meson exchange model for es- 
timating the dissociation of charmonia in the hadron gas 
and a quasifrcc particle approximation for charm quark 
inside the J/tj} in the quark-gluon plasma. Binding ener- 
gies of charmonia become very small at high temperature 
so that the charm and anti-charm quarks inside the char- 
monia indeed can be approximated by a quasi free par- 
ticles. With this idea, GR approximate the dissociation 
cross section of charmonia with the elastic cross section 
of charm or anti-charm quark, once the energy transfer 
is larger than the binding energy of charmonium. 

However, as we show in the appendix, such contribu- 
tion constitutes only the leading monopolc contribution, 
whose contribution from the charm and anticharm quarks 
cancel as they have opposite charge. A consistent calcu- 
lation shows that the dissociation cross section is of the 
dipole type and only sensitive to the size of the wave 
function (r 2 ) 1//2 as is expected of a system composed of 
a quark and an antiquark system with opposite charge. 



THERMAL DECAY OF CHARMONIA 



B. Thermal width in the hadron gas 



The thermal decay width of charmonia at temperature 
T can be calculated using the following factorization for- 
mula (H: 



The dissociation cross section of charmonium by 
hadron i can be calculated by the following factorization 
formula: 
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°f ss (s) = Yl [ dx nij (x,Q 2 ), (18) 

3 

where riij(x, Q 2 ) is the distribution function of parton j 
in hadron i. x is the longitudinal momentum fraction 
of parton j in hadron i, which is a value between and 
1. Q is the renormalization scale of the parton distri- 
bution function. Suppose that a gluon is emitted from 
a light quark in a hadron. If transverse momentum of 
the gluon is smaller than Q, parton distribution function 
absorbs this splitting process. But if the transverse mo- 
mentum is larger than Q, this splitting process has to 
be calculated in <7j(xs,Q 2 ) of Eq. (fl8]) . In other word, 
dissociation cross section of charmonia by hadron i, the 
left hand side of Eq. (TT8")) . is composed of nonperturba- 
tive part, which is the parton distribution function, and 
perturbative part, which is the dissociation cross section 
of charmonia by parton j. That is why Eq. (fT8|) is called 
factorization formula. 

The factorization formula also shows how collinear di- 
vergence is removed. If a massless gluon is emitted in the 
same or opposite direction from the original light quark 
or a massless gluon, then the propagator will have a pole 
and induce a collinear divergence. In this case, trans- 
verse momentum of the emitted gluon is smaller than Q 
and thus the divergence can be absorbed by the parton 
distribution function. The elementary cross section Uj as 
well as the parton distribution function mj depends on 
the scale Q 2 but the final result, af lss does not depend 
on the scale in principle. There are several parton distri- 
bution functions available. However the minimum scales 
Q above which they are defined, are all larger than the 
binding energies of charmonia, which is the renormaliza- 
tion scale for aj(xs,Q 2 ). For the case of heavy quark 
scattering, the separation scale Q can be taken to be the 
transferred energy or the mass of heavy quark. But in 
the case of heavy quarkonium, the separation scale is the 
binding energy of heavy quarkonium. 

We use the Bcthc-Salpetcr amplitude to describe the 
bound state of quarkonium in which the generated po- 
tential becomes Coulombic in the heavy quark limit; in 
principle, a consistent counting and renormalization is 
possible in such limit. Although the phenomenological 
heavy quark potential for strong force is composed of the 
Coulombic potential part and the linearly rising string 
potential energy [3l|, in the heavy quark limit, the dis- 
tance between quark and antiquark is very small and the 
effect of the linearly rising string part of the potential 
is not felt by the quarks. In such a limit, the binding 
energy of J/ip, which is IS state, can be estimated using 
the mass difference to its first excited 2S state, which will 
also be Coulombic. Although the charm quark mass is 
not so large, we can still approximate the states to be 
Coulombic and obtained the binding energy of J/ip to be 
780 MeV [H, H3] ■ Unfortunately, no parton distribution 
functions are defined in such a low Q. 



For our purpose, we just take the MRSSI parton dis- 
tribution function [HJ of the pion at V5 GcV and use 
it for all hadrons as the pion is the most abundant par- 
ticle in the hadron gas stage. The minimum scale Q of 
MRSSI is still much larger than the binding energy of 
J/ip- In fact, due to this reason, it was found that com- 
bining the parton distribution function at larger Q value 
than the scale of the scattering cross section in Eq. [TS] 
we found inconsistent results for the charmonium states 
where the cross section becomes negative at some incom- 
ing energies [3(3 . Therefore, we will first apply the ap- 
proach to calculate the dissociation cross section for the 
bottom system, perturbative QCD approach always gives 
positive cross section[3(| and then extrapolate the result 
to the J/ip case by the following dipole approximation, 

^U/i)=(^)Vf(v^). (19) 

Here, i?j/^,,i?x are the radii of J/ip and T respectively. 

In Coulombic potential, the radius R = l/^me , 
where m is mass of constituent heavy quark of quarko- 
nium, 1.94 GeV for charm, 5.1 GeV for bottom, and e 
is the binding energy of quarkonium. The difference in 
the renormalization scales of the couplings appearing in 
the J/ip and T cross section appearing in Eq. (|19j) is ne- 
glected. Two different energy scales y/s, y/s are assumed 
to have the relation yfp = (m c /mb)y/s, where m c , mb are 
masses of constituent charm and bottom quarks. 

The upper figure of Fig. [6] shows the dissociation cross 
section of J / ip by the tt obtained using the method above. 
The present cross section is still smaller than that esti- 
mated using the meson-exchange model [3j| HH or the 
quark-exchange model [37j • 

The cross section obtained above is now put to Eq. 
(fT6|) to get the thermal width of J/ijj by tt. For more 
simplicity, the thermal width in hadron gas is simplified 
as following: 

T(T) = Aj ^v rel (k)n^(k,T)a diss (k,T), (20) 

where v re i is relative velocity between charmonium and 
pion in rest frame of charmonium, and a dtss is dissocia- 
tion cross section of charmonium by the pion. A is the 
ratio of number density integrated in momentum space 
of all hadrons, whose masses are less than 1.5 GeV for 
mesons and 2.0 GeV for baryons, and that of pion. The 
integrated number density of pion is about 0.16 /m~ 3 
at temperature of chemical freeze-out, and that of all 
hadrons is about 0.43 fm~ 3 . 

In order to calculate the dissociation cross section of 
charmonia in perturbative QCD approach, we have to 
know their wavefunctions [29l . |30| . For J/tp in hadron 
gas , wavefunction of IS Coulombic bound state is used 
[30(. For excited charmonia in hadron gas, it is assumed 
that the cross section is proportional to squared radius. 
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FIG. 6: Upper figure: the dissociation cross section of J/ip 
by 7r. Lower figure: the thermal widths of J/i/>(solid line), 
\c (dotted line) and i// (dashed line) in the hadron gas and in 
the quark-gluon plasma obtained in perturbative QCD ap- 
proach. 



According to the screened Cornell potential [3l|, vac- 
uum screening masses for J/tp, Xc an d ip' are 0.18, 0.18, 
and 0.26 GeV, from which the radii are found to be 
0.563, 0.778, and 1.504 fm respectively. As a result, 
the dissociation cross section of Xc (VO is 1-9 (7-1) times 
larger than that of J/tp in the hadron gas. 
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C. Thermal width in the QGP 

The thermal width depends on the charmonium wave 
function. Here we use the Cornell potential with temper- 
ature dependent screening mass. The form of screened 
Cornell potential is given as 

V(r,T) = ^(l-e-> lr j -^e-^ r , (21) 

where a is 0.192 GeV 2 and a = 0.471. The screening 
mass fi in quark-gluon plasma depends on temperature. 
In the limit that fi goes to zero, the screened potential 
becomes the normal Cornell potential. Fig. [7] shows the 
wavefunctions of charmonia in momentum space at vari- 
ous screening mass. We can see that the wavefunctions of 



FIG. 7: Wavefunctions of J/ip (upper), Xc (middle), and ip' 
(lower) in momentum space, at screening masses, 289 MeV 
(solid line), 306 MeV (dotted line), and 323 MeV(dashed 
line) , and 340 MeV (dot-dashed line) . 

Xc & n d ip' are sensitive to the change of screening mass, 
while that of J/tp is not and is more localized at the 
origin in momentum space. 

For the temperature dependence of the screening mass, 
we simply use the leading order pQCD result, /i = 
y/(N c /3) + {NfJ&)gT. The coupling constant g is 1.47 in 
the case of quark-gluon plasma, which is obtained from 
the best fit for Raa data in our formalism. The details 
for g in quark-gluon plasma will be explained later. In the 
case of quark-gluon plasma, the dissociation cross section 
does not have divergence, because the effective thermal 
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masses of partons play the role of cutoff in momentum 
space; in a sense the quasi particles are the asymptotic 
states of the QGP that can perturbatively scatter with 
the charmonium states. As long as the temperature scale 
is smaller than the separation scale, such description can 
be thought of as an effective factorization formula, where 
all the soft physics is included in the quasi particles. 

The lower figure of Fig. [B] shows the thermal widths of 
charmonia both in the hadron gas and in the quark-gluon 
plasma, calculated using our NLO formula. We see that 
the thermal widths of charmonia suddenly increase by 
about an order of magnitude at the critical temperature. 
There are several reasons for such a sudden increase of the 
thermal widths. One of them is the difference in kinetic 
energies of partons. The kinetic energy of a parton has 
to be larger than the binding energy of the charmonium 
states for the break up to take place. Such energy are 
readily available for massive quarks and gluons at finite 
temperature. However, for pions, the temperatures have 
to be higher. Moreover, there is a large increase in the 
degree of freedom as phase transition occurs. 

Another important factor for survival rate of J/ip is 
the feed-down contribution from its excited states such 
as the Xci ip' ■ Some J/tp is produced directly but some of 
them is produced through the decay of its excited states. 
It is not clear yet what fractions of J/tp are produced 
through feed-down. Here it is assumed that 25% and 8% 
of J/tp are fed down from Xc and tp' respectively, which 
is the world average [IH . Then the survival rate of J/tp 
from thermal decay is 

S th (b,s) = 0.67S^(b,s) + 0.25S%(b,s) 

+0.08S%(b,s), (22) 

where 

S{ h (b, s) = cxpj - jf " P(r')dr'\ (s e W) 

= (s £H J ). 

The superscript j in the lower equation is species of 
charmonia, and r c /, the upper limit of the integration, 
is proper time for the chemical freeze-out. H° is the hot 
region where initial temperature is over the dissociation 
temperature of each charmonia. 

VI. THE REGENERATION OF CHARMONIA 

As discussed before, results from hydrodynamics sug- 
gest that hot nuclear matter created in relativistic heavy 
ion collision is thermalized very fast. Moreover particle 
yields and their ratios from relativistic heavy ion collision 
are well reproduced using the grand canonical ensemble. 
If the number of charm and anti-charm quarks are also 
thermalized, the total number of charm pairs produced 
N^ B (t) from the collision of nucleus A and B would sat- 
isfy the following equation : 



N^{r) = |in (r) + n,(r)j V». (23) 

Here, n a and are the number density of open and 
hidden charms respectively, and V is the volume of hot 
nuclear matter; all are a function of the proper time r. 

In general, if creation or annihilation of charm pairs is 
considered during the expansion of hot nuclear matter, 
the number of charm quark pairs will be a function of 
proper time r. Moreover, if the lifetime of hot nuclear 
matter is sufficiently longer than the thermalization time, 
the number of charm pairs will become thermalized and 
satisfy Eq. (|2"3")k thermalization will proceed through the 
processes such as c + c o q + q, c + c •<->• g + g. How- 
ever, estimates based on perturbative QCD suggests that 
the lifetime of hot nuclear matter is not long enough for 
the number of charm pairs to be thermlized, because the 
cross section for creation or annihilation of charm pairs 
is small with respect to the life time of the fireball. 

Therefore, taking into account the off-equilibrium ef- 
fects for the heavy quarks, the number of charm pairs 
given in Eq. ([2U)) will be modified as follows, 

N* B (t) = |^o(r) + 7 V(r) jv(r), (24) 

where 7 is called the fugacity Q . If the number of charms 
is more (less) then what it would be if the charm quarks 
are thermalized, fugacity is larger (smaller) than unity. 
In the central collision of RHIC the initial temperature of 
created hot nuclear matter is very high and the fugacity 
at t = To is smaller than 1. As the system cools down, 
the fugacity increases. This fact can in fact be simply 
understood qualitatively as follows; assuming isentropic 
expansion, V(t) cx I/ul where ul is the density of light 
partons, substituting this into Eq. (f24)) gives N^ 5 b {t) oc 

|l'T 2 tr^+^ 2! ^r^|' since ]h ^T ~ exp(-m c /T), where 

m c is the mass of charm quark, and assuming the number 
of charm quarks do not change during the expansion, 
we can conclude that the fugacity should increase as the 
temperature decreases during the expansion. The final 
value of fugacity at T c as a function of A^^t is shown in 
the lower figure of Fig. [8] 

In Eq. (|24|) the number density of open charm is multi- 
plied by 7 and hidden charm by 7 2 , because open charm 
has one charm or anti-charm quark while hidden charm 
has both charm and anti-charm quarks. It may be as- 
sumed that the number of charm pairs does not change 
during the expansion of hot nuclear matter, because the 
creation and annihilation cross section of the pair is very 
small. Then the left hand side of Eq. (|24|) dose not 
change and its initial value is calculable with Glauber 
model as following: 

N* B (b) = v» N AB J d 2 sT A {s)T B {b-s) 

= ag N ABT AB (b), (25) 
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where u^ e N is the cross section for cc pair production in 
p+p collision, which is 63.7 [ib per rapidity near midra- 
pidity at y/s = 200 GeV in perturbative QCD (H|33|. 

Another correction is about system size. So far grand- 
canonical ensemble is used for number densities of parti- 
cles. Canonical ensemble, however, is to be used techni- 
cally, because quantum numbers of the system arc to be 
conserved. Quantum numbers of the system are fixed in 
canonical ensemble, while they are not in grandcanonical 
ensemble. The function of probability density with re- 
spect to quantum number is a delta function in canonical 
ensemble, while it is a gaussian type distribution in the 
grandcanonical ensemble. As the system size is large, the 
width of function of probability density in grandcanon- 
ical ensemble is narrow. In the thermodynamical limit, 
where the system size is infinite, the function becomes 
like a delta function. Therefore, grandcanonical ensem- 
ble and canonical ensemble are equal in thermodynamical 
limit. The equality seems to be true in central collisions 
of heavy ions, because the number of produced charm 
or anticharm are not small. However it is questionable 
for peripheral collisions. The quantum number in grand- 
canonical ensemble can be converted to that in canonical 
ensemble by simply multiplying modification factors as 
follows; 



N, 



AB 



1 ^ h( lno (T)V{T)) 

2 in ° [T) Iohn (r)V{T)) 

+ 7 2 n h (r)W(r) 



(26) 



where I and I\ are modified Bcssel functions and their 
ratio that depends on the number of charm quarks is 
called the canonical ensemble correction 

The upper figure of Fig. [5] is the canonical ensem- 
ble corrections at hadronization. In central collisions, 
the number of produced charm quarks are not small so 
that the correction is almost unity. However the cor- 
rection rapidly decreases in peripheral collisions because 
the number of produced charm quarks becomes small. 
In other words, the difference between grandcanonical 
ensemble and canonical ensemble becomes serious. 

Fugacity 7 is calculated from Eq. (|25)l and ((26)) . The 
lower figure of Fig. [5] shows the fugacity at hadroniza- 
tion. The Fugacity is about 6 at central collision and 
remains constant until N part decreases to about 100, and 
then rapidly increases. The rapid decrease of canonical 
ensemble correction causes fugacity to rapidly increase at 
small N part . 

So far it is assumed that charm quarks are thermal- 
ized kinetically. In other words, it is assumed that the 
kinetic energy of the charm quarks follow thermal distri- 
bution. This thermalization is attained through elastic 
scattering of charm quarks with the surrounding partons. 
However the elastic, as well as the inelastic, cross sections 
of the heavy quarks are smaller than those of the light 
quarks. The formation of charmonium states through 
recombination are more probable when the charm and 
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FIG. 8: Canonical suppression factor(upper) and fugac- 
ity (lower) with respect to number of participants 



anti-charm quarks are randomly distributed in momen- 
tum space. Therefore, to estimate the regeneration ef- 
fect, it is important to know how much charm quarks 
and anti-charm quarks are thermalized kinetically at the 
hadronization point. 

To estimate this effect, one defines the relaxation factor 
as follows, (hi : 



R = 1 — cxp 



T QGP 



(27) 



' <::</ 



where tqqp is the proper time the quark-gluon plasma 
phase terminates. Relaxation time, r eq , is defined as fol- 
lows: 



*CT) = i/(E/ 



dk 3 
(M 



(fc)ni(fc, , (28) 



where n,; is the number density of parton i in the quark- 
gluon plasma, and crj is the elastic cross section of charm 
or anti-charm quark by parton i in the charm or anti- 
charm rest frame. The cross section multiplied by num- 
ber density of partons is the mean free path of charm or 
anti-charm quark. Therefore r eq is the average elastic col- 
lision time of charm or anti-charm quarks with the light 
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partons. The clastic cross section is calculated to leading 
order in perturbative QCD. The effective thermal masses 
of partons extracted from lattice by Levi and Heinz are 
used [IH, H^. The relaxation factor R is about 0.25 at 
central collision and decreases as the collision becomes 
peripheral. 

The number of regenerated J/ip can now be written as 

KisW = l 2 {nj/^ H + Br( Xc ^J/^)n Xc S%_ H 

+ Br{i>' -> J/^)n^S^_ H JkR, (29) 

where Br(x c — > J/ 'VO and Br(tp' — > J/ip) are branching 
ratios from \c to J/ip and from ip' to J/ip respectively. 
Canonical ensemble correction is not multiplied, because 
charmonia are hidden charms. S^I^_ H , Sfu_ H and Sf h _ H 
are respectively the survival rates of J/ip, Xc and ip' from 
thermal decay in hadron gas defined as follows: 

Sk-H = cxpj - jT" r l (r')dr'|. (30) 

Here, tr is the proper time the hadronization is com- 
pleted. 

VII. RESULTS 

J/ip is a massive particle produced through hard colli- 
sion of two nucleons from each colliding nucleus. If there 
is no nuclear modification, the number of produced J/tp 
in a heavy ion collision will be proportional to the num- 
ber of binary collisions as follows, 
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FIG. 9: The upper figure compares Raa from the experiment 
[l^ ] to that calculated from Eq. <|33p . The dashed line is 
the nuclear absorption and thermal decay, the dotted line the 
regeneration effect and the solid line their sum. As for the 
lower figure, the experimental data and the solid line is the 
same as in the upper figure, but the dashed line shows the 
result when the initial suppression of J/ip in the hot region is 
not taken into account. 



Nf}^(b) = a^T AA (b) 



(31) 



Here, is the cross section for J/ip production in p+p 
collision, which is recently measured to be 0.774 ub per 
rapidity near midrapidity at ^fl = 200 GeV [H, El| . Taa 
is the thickness function defined in Eq. ([T|). 

The nuclear modification factor Raa is defined as fol- 
lows, 



Raa{V) 



Nf/ A 4b) 



aJ^T A A(b) 



(32) 



The numerator in the right hand side of the above 
equation is the actual number of J / ip produced in A+ A 
collision, while the denominator is the expected num- 
ber of J/ip production in the same collision. However, 
nuclear absorption and thermal decay of charmonia sup- 
press the actual number of produced J/ip so that these 
effects supress Raa- On the contrary, regeneration of 
charmonia enhances the number so that it raises Raa- 



Considering these effects altogether, the nuclear modifi- 
cation factor is from Eq. g]), and (1521 



R AA {b) 



d 2 sS„ 



2 (b,s)S t h(b,s) 



N AA 
1S reg J/4, 



(b) 



a™A*T AA {b) 



(33) 



The first term of the right hand side represents nuclear 
absorption and thermal decay of J/ip, which includes 
the convolution in the transverse plane that takes into 
account the initial suppression if the local temperature 
is larger than the charmonium dissociation temperature, 
and the second term regeneration of J /ip. 



A. results for Au+Au 

The upper figure of Fig. [3] shows that Eq. (|3"3")l , with 
a best fit coupling q = 1.47, reproduces the experimen- 
tal data of Raa [131 weu - The dashed line is S nuc S t h, 
the dotted line R^a an d the solid line their sum. The 
coupling constant g is an important parameter in the 
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result. It determines the thermal widths of charmonia, 
their dissociation temperatures and the relaxation fac- 
tor R. When the coupling constant increases, the ther- 
mal widths of charmonia become wider, the dissociation 
temperatures of charmonia become lower and the relax- 
ation factor R larger, and vice versa. While perturbative 
QCD formulas are used, it is known that the quark-gluon 
plasma from T c up to about 2 T c is nonperturbative. In 
this sense, the coupling constant g is kind of a free pa- 
rameter fitted to reproduce the experimental data. 

In our approach the value of the coupling g also deter- 
mines the position of a kink on Raa curve. At present, 
the experimental uncertainties in Raa is still too large 
to extract any detailed structure as a function in -/V part . 
Nevertheless, taking the central experimental data of 
Raa seriously, there seems to be two drops. The first 
one appears at small number of participants, and the 
second one between 170 and 210 N part . 

There are two important thermal effect on charmonia 
production: one is the initial suppression which occurs 
when the local initial temperature of the hot region is 
higher than the dissociation temperature, the other is 
thermal decay of charmonia through interaction with sur- 
rounding particles. Two effects occur at different temper- 
atures, which are well distinguished in the lower figure of 
Fig. |U From the critical temperature to the dissociation 
temperature, which are respectively identified by the sud- 
den increase and the divergence of the thermal width, the 
thermal width of charmonium increases monotonically; in 
the lower figure of Fig. HI the divergence is represented by 
the termination point. For the initial suppression effect, 
the survival rate of charmonium is 0, while for the effect 
coming from thermal decay, it depends on the tempera- 
ture and lifetime of the QGP relevant for thermal decay. 

It should be noted that the survival rate will change 
appreciably only when the thermal width becomes large, 
because the lifetime of quark-gluon plasma in heavy ion 
collision is not so long. In fact, as can be seen from 
Fig. [5j the temperature region near T c is more impor- 
tant, which last almost 2 to 3 fm/c. This means that the 
thermal width has to be of the order 100 MeV to have an 
appreciable suppression effect. Therefore, there will be a 
critical change in the survival rate depending on whether 
or not the initial temperature is larger than the disso- 
ciation temperature; because that condition corresponds 
to having either an infinite or finite thermal width. This 
critical change in survival rate will appear as a sudden 
change of Raa along the axis of number of participants, 
because the initial temperature rises as the number of 
participants increases. Therefore, if at some future point , 
the error bars in the experimental data are reduced and 
one indeed finds the sudden kink, it would strongly im- 
ply that the J/ip dissociation takes place at the corre- 
sponding point. Moreover, an initial kink could appear 
when the initial temperature is higher than the dissocia- 
tion temperatures of excited charmonia, such as \c and 
ip'. It should be kept in mind however that when the 
thermal decay width is in the order of few hundred MeV 



just above the deconfincment, the sudden kink could also 
have occurred due to the thermal decay. 

Again, taking the central values of Raa seriously, the 
second kink can be translated to a dissociation tempera- 
ture of 386 MeV; this can be obtained from Fig. [2] and 
the Glauber model given in Eq. ((T|) that links A part to 
the impact parameter. Now, solving Schrodingcr equa- 
tion with the screened Cornell potential, one finds that 
the dissociation occurs when the screening mass is about 
fx = 695 MeV. Since the screening mass is simply 
\x = y/(N c /3) + (Nf/6)gT in our calculation, the cou- 
pling constant g has to be fixed at 1.47 to reproduce the 
screening mass at the dissociation temperature 386 MeV. 
This coupling relates to a = 0.172 which is not so differ- 
ent from the lattice calculation for the asymptotic value 
of the coupling reached at higher temperature [42j and to 
the coupling given in Eq. (fTTI) . 

With the same value of g, one finds the dissociation 
temperatures of Xc and ip' to be 199 MeV and 185 MeV 
respectively. The initial drop of the solid line in both 
the upper and lower figures of Fig. [5] is coming from the 
initial suppression of both \c and ip'. The dashed line of 
the lower figure of Fig. [9] corresponds to the case where 
the initial suppression of J /ip is not taken into account. 
We can see how much Raa is improved by considering 
the effects of the initial hot region. 



B. understanding the results 

To gain some more insights into what causes the char- 
acteristic features of the solid line in Fig. O Let us 
artificially change some parameters and investigate the 
changes in the prediction. The results are shown in 
Fig. [lOl The solid line is the best fit value using 
the same method as described above but obtained with 
g = 1.47 and with a constant thermal width of 10 MeV. 
The dashed and dot-dashed curves are when the ther- 
mal width is changed to 100 MeV and zero respectively. 
One sees that the slope of the curve critically depends on 
the thermal width. While the increase in the width sup- 
presses Raa, one should also note that an increase in the 
width also implies a larger coupling to the medium and 
hence an increase in the relaxation factor. This will en- 
hance the regeneration effects and thus compensate the 
suppression. 

In the dotted curve, the mass of J/tp is artificially 
changed by -100 MeV. This is to see if the sudden mass 
shift of J/tp recently advocated by Morita and one of 
us[5| will have any effect on Raa- As one can see, due to 
the enhancement of the statistical factor in regeneration 
effect, the curve is pushed up by a nontrivial amount. 
Therefore, the effects of a sudden decrease in the mass 
and a sudden increase in the width at the phase transition 
point, as advocated in Rcf. 0, |43|, might just compen- 
sate each other in Raa, making it difficult to identify 
such effects from a measurement of Raa alone. 
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FIG. 10: Raa when the thermal width and mass are artifi- 
cially changed. The different curves are explained in the text. 



C. results for Cu+Cu 

The same method is applied to Cu+Cu collision at 
\/Snn = 200 GeV. The only difference from the previ- 
ous analysis is the geometry of the colliding nuclei. r a 
and C for Cu are set at 4.163 fm and 0.606 fm respec- 
tively in Eq. ([3]) [2l[. The upper figure of Fig. [TT] shows 
Raa of J ftp in Cu+Cu collision at ^/Snn = 200 GeV. 
Transverse acceleration of hot nuclear matter in Cu+Cu 
collision will be less than that in Au+Au collision, so 
that two different transverse accelerations, 0.1 c 2 /fm and 
0.08 c 2 /fm, are used in the calculation for Raa- How- 
ever there is almost no difference in Raa so that the two 
curves almost overlap in Fig. [TTJ The calculation for the 
Cu+Cu collision slightly overestimates the experimental 
data at central collisions while it underestimates at pe- 
ripheral collisions. In Cu+Cu collision at \/Snn = 200 
GeV, the maximum temperature is lower than the dis- 
sociation temperature of J/?p even at the most central 
collision, hence no kink appears in Raa- 
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FIG. 11: The upper figure is Raa of J/tp in Cu+Cu collision 
at y/SNN = 200 GeV. The solid line for a± = 0.1 c 2 / fm and 
the dotted line for a± = 0.08 c 2 / fm are overlapped. The lower 
figure is Raa of J/tp in Pb+Pb collision at ^/snn = 5.5 TeV. 
The dashed line is Raa after nuclear absorption and thermal 
decay, the dot-dashed line is Raa from regeneration effect and 
the solid line is their sum. 



D. results for LHC 

We can attempt to predict Raa of J/ ip in Pb+Pb col- 
lision at LHC as well. Because the collision energy is 
much higher than at RHIC, parameters will have to be 
changed from those used before. We change only some of 
them. r a and C for Au are set at 6.62 fm and 0.546 fm 
respectively [2l|. 6.4 /ib and 639 /ib are use for the cross 
sections for J/tp production and for cc pair production 
per rapidity in p+p collision (T^]. The constant 30.3 for 
produced entropy in Eq. §5§ is replaced by 55.7 [HI]. The 
variable x in Eq. |(5J) will be larger than 0.11, because x 
increases as collision energy increases. However we use 
the same value for lack of information. All chemical po- 
tentials are set at zero, because midrapidity region will 
almost be baryon-free at LHC energy. All other param- 
eters including transverse acceleration and its terminal 
velocity of hot nuclear matter are set at the same values 
as at RHIC for simplicity and for lack of knowledge. 

The lower figure of Fig. [TTJ shows the prediction for 
Raa at LHC energy. In contrast to the case at RHIC en- 
ergy, the survival rate from thermal decay is very low and 
the regeneration effect is dominant except in the periph- 
eral collision. This is so because the lifetime of quark- 
gluon plasma is much longer than at RHIC. According 
to our rough estimation, the lifetime at the most cen- 
tral collision is twice that of the RHIC case. The long 
lifetime of quark-gluon plasma suppresses the survival 
rate of charmonia through thermal decay and enhances 
the relaxation factor R of charm and anti-charm quarks. 
Abundance of charm quarks is another reason for the 
dominance of the regeneration effect. 

The same method seems to be applicable to lower col- 
lision energies such as SPS. In this case, however, the 
baryon chemical potential is too high to make use of ef- 
fective thermal masses by Levai and Heniz, which are 
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extracted from the lattice QCD data with zero baryon 
chemical potential. 



VIII. CONCLUSION 

We introduce a generalized two-component model to 
calculate the Raa of J /if) near midrapidity in Au+Au 
collision at y/s = 200 GcV. In the model, suppression is 
caused by the nuclear absorption of primordial charmonia 
and by the thermal decay of charmonia in the quark- 
gluon plasma and in the hadron gas. Enhancement is 
caused by the regeneration from the QGP. One of our 
emphases is the use of a consistent perturbativc QCD 
approach, up to next to the leading order, to estimate 
the thermal decay widths of charmonia both in the quark- 
gluon plasma and in the hadron gas. While the coupling 
is considered to be a free parameter to be determined by 
experiment, it is also related to the screening mass and 
subsequently dissociation temperatures of charmonia in 
screened Cornell potential model. 

Another emphasis of our approach is the initial tem- 
perature profile which leads to initial suppression of char- 
monia if the local temperature is higher than the disso- 
ciation temperature. When the initial temperature pro- 
files, and thus the initial hot regions are considered, the 
calculated Raa seems to show two sudden drops; one 
at very small number of participants and the other at 
around 170 ~ 210. The first sudden drop of Raa at very 
small number of participants is caused when the maxi- 
mum temperature of hot nuclear matter reaches above 
the dissociation temperatures of excited charmonia and 
the second sudden drop when it reaches the dissociation 
temperature of J/ip. On the other hand, if the thermal 
width of J / ip is large compared to the inverse time scale 
of the QGP, the second sudden drop of Raa disappears. 
Moreover, some effects can cancel other effects. For ex- 
ample, a increased thermal width of J/ip will reduce the 
Raa, while a decreased mass of J/ip will enhance it. 

The present experimental data from RHIC are still 
large to discriminate all the different and competing 
mechanisms discussed here. However, the central val- 
ues seem to be consistent with the two drops. In this 
sense, the maximum temperature of SPS seems to be 
lower than the melting temperature oiJhp, because the 
second drop is not seen in Raa there [44| . Similar con- 
clusion could also be made from the data for the Cu+Cu 
collision at yfs = 200 GcV. However refined data for the 
Au+Au collision at y/s = 200 GeV are required to make 
a definite conclusion. Moreover fitting a single data can 
not disentangle all the different suppression mechanisms 
separately. In that sense, the Raa of J/ip at LHC en- 
ergy should provide a very interesting information. As 
shown in our rough estimate, Raa at LHC is very dif- 
ferent from that at RHIC, in that it is dominated by the 
regenerated J/ip from the QGP. Moreover, within the 
thermal decay width given in the paper, the life time of 
the plasma seems to be much longer than the inverse of 




M k2 



FIG. 12: one diagram of quark- induced J /if) dissociation at 
next to the leading order in perturbative QCD, and its de- 
composition 

the thermal decay with. Hence, initially formed J/tp will 
be mostly suppressed through thermal decay. Therefore, 
the upcoming LHC experiment will provide a totally new 
information about J/ip in relativistic heavy ion collision 
and thus the nature of QGP and its formation. 
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Appendix 

The upper figure of Fig. [TJ] is a diagram for the J/ip 
dissociation by a light quark at next to the leading order 
in perturbative QCD. The double line represents exter- 
nal line of J /ip and the adjacent small circle the Bethe- 
Salpctcr amplitude. The Bcthc-Salpctcr vertex repre- 
sents bound state of charm quark and anti-charm quark 
and has the following form in the heavy quark limit and 
in rest frame of J/ip: 



where [i is the spin index of J/ip and e its binding energy. 
p is the relative three momentum of p\ and P2 , p = \p[ — 
P2I/2, and ip(p) the wavefunction of J/tp. N c is number 
of color. The charm quark propagator between Bcthe- 
Salpctcr amplitude and gluon vertex is reexpressed as 
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k 2 



m c _ . J2s uS ( k ) uS ( k ) 
m 2 k 2 — m 2 



(35) 



If the binding energy of J/ij; is small, the momentum 
k of charm quark propagator will be slightly off-shell. By 
using Eq. (f3"5j) . the invariant amplitude for Fig. [T2"lis 



(ki - k 2 ) 2 - m; 



■u(pi)j u u s (k) 



xu s (k)T l _ L (k,p 2 )v(p 2 )- 



Y,MQFU s {k)Y^k,~p 2 )v{p 2 ) 



k 2 



(36) 



where Mq F is the invariant amplitude for c + <7 — > c + q 
with spin s of incoming charm quark in the approxima- 
tion of quasifree particle The lower figure of Fig. [T21 
shows the diagrammatic decomposition of Eq. (f36|) . In 
heavy quark limit, the condition of energy conservation, 



or 



h o_ k o_ f _ \fi\ 2 , \P*\ 

™1 "<2 ° — 



2m2 2m 2 ' 

gives the following order counting 

fc° ~ fc 2 ° ~ 0(mg A ) ~ |p 2 | ~ 0(m 3 2 ), (37) 

because the binding energy e Q is of mg 4 order [3(| ■ Un- 
der this order counting, the denominator of charm quark 
propagator can be approximated as follows, 

k 2 — m 2 c = {pi — k\ + k 2 ) 2 — m 2 

rj -2pi • (fci - fc 2 ) « -2m c (fc? - 

9 / JpiTJpII 2 



JpT 

-2m c e H - 

mi 



(38) 



where 



p = (pi - P2)/2 = pi - (fci - fc 2 )/2 « pi, 

p = (pi - pl)/2 = -pi + (fci - fe)/2 ~ -pi 



are used for the last approximation in Eq. ([55)1 . The 
denominator of charm quark propagator indicates how 
off-shell the propagator is. We can see the off-shellness is 
canceled by the same factor in Bethe-Salpetcr amplitude. 
Therefore the invariant amplitude for J/ip dissociation is 



= —2i m c 



n J/j) 



4>{p)Y. M QF 



i + y . i - y 

xw s (fc)^-7,5^^-«(P2), (39) 

and its averaged square with respect to the spin and color 
of incoming particles is 



|Mr = 



1W 2 



\i>(p)\ 2 \M QF f 



= 2m J/ ^{p)\ 2 \M QF \ 2 . (40) 

The phase space of final state is separated as follows, 

d 3 p\d 3 p 2 d 3 k 2 
(27r) 3 2 J Bi(27r) 3 2 J B 2 (27r) 3 2fcO 

x {2ir) 4 5 4 (p 1 +p 2 + k 2 -q-k 1 ) 



d 3 p 2 d A k 



{2n) 3 2E 2 

d 3 pxd 3 k 2 



(2TT) 3 2E 1 (2Tr) 3 2k° 



S\k+p 2 -q) 

(2Tr) 4 S 4 ( Pl +k 2 -k-k 1 ). 



(41) 



Finally, the dissociation cross section of J/ip from the 
diagram is 



diss 



1 



d 3 p\d 3 p 2 d 3 k 2 



Amj/^kl J (2n) 3 2E 1 {2n) 3 2E 2 (2n) 3 2k 
x (2ir) i 6 4 (p 1 + P2 + k 2 -q- fci)|M| 2 



d 3 p 2 d A k 4/ 
{2n) 3 2E 2 



1 



-8\k+p 2 -q)\4,(j>)Y 
d 3 pid 3 k 2 



2fc° J (27r) 3 2E 1 (2n) 3 2k° 
x<5 4 (pi+fc 2 -fc-fci)|Mo7| 2 . (42) 

In heavy quark limit, 

E 2 mm c , k = q - p 2 tn (m c , 0), (43) 

the dissociation cross section is simplified as follows, 



d 3 p 



|V'(P)|V 



2 „QF 



(44) 



where p 2 is approximated to —p. <jQ f is the cross sec- 
tion of charm quark, which is used for dissociation cross 
section of J/ip in the approximation of quasifree parti- 
cle. Suppose that the cross section is independent of 
p, which enters in <jQ f through the momentum carried 
by the quasifree charm quark, then the momentum inte- 
gral with the square of the wave function appearing in 
Eq. (|44[) just factors out as the wave function normaliza- 
tion, which is unity, and we will have 



a diss = a® F . 



(45) 
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FIG. 13: other diagrams of quark-induced J/tp dissociation 
at next to the leading order in perturbative QCD 



so they are treated equally. It gives the following result 



IMI 2 = - 



g i m 2 c m J/4: 



dip(p) 



dp 

1 (*0)2_, 

2 



(fc 2 °) 5 



2fci • fc 2 



(46) 



In this limit, the cross sections will have to be just a 
constant independent of the incoming energy. 

However, Fig. [12] is not the only diagram for quark- 
induced J /tp dissociation. The left figure of Fig. [T3]is an- 
other diagram for the same dissociation. In the quasifree 
particle approach, this diagram corresponds to the pro- 
cess c + q — ?> c + q while the diagram of Fig. [12] to the 
process c + q — > c + q. Moreover, the invariant ampli- 
tude for each process is squared and then summed to get 
dissociation cross section. As a result, cross section be- 
comes double the single contribution. On the other hand, 
in our perturbative QCD approach, the amplitude of the 
two diagrams are summed before being squared. In this 
case, the invariant amplitude for the left figure of Fig. 
[T51 cancels that of Fig. [T2] exactly so that the next to 
the leading order in the counting scheme of Eq. (|37[) has 
to be considered. In the order counting, the next to the 
leading order of Fig. [12] the left figure of Fig. [13] and 
the leading order of the right figure of Q2] are all the same 



It is worthy to notice that Eq. (|46l) is a function of 
dip(p)/dp, the dipole of wavefunction, while Eq. (|l4"]l is 
a function of ip{p)j the monopole of wavefunction. The 
result with the dipole form is natural result in the sense 
that J/\jj is a dipole of color charge. 

The comparison between two models in the gluon- 
induced dissociation cross section is similar. The only 
difference is that perturbative QCD approach has dia- 
grams where one gluon is absorbed at a heavy quark 
or antiquark line and the other gluon is emitted at the 
other heavy antiquark or quark line as shown in dia- 
grams (3) and (4) of Fig. 6 in [13, but the approach 
of quasifree particle doesn't have such diagrams. In ad- 
dition to the differences discussed above, comparing the 
quasifree model with ours, one finds that the perturba- 
tive QCD approach has contributions from the diagrams 
where gluon is emitted or absorbed at the gluon line 
exchanged between the heavy quar k and antiquark as 
shown in (9) - (16) of Fig. 6 in ip. 
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